{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');\n",
       "@import url('http://fonts.googleapis.com/css?family=Vollkorn');\n",
       "@import url('http://fonts.googleapis.com/css?family=Arimo');\n",
       "@import url('http://fonts.googleapis.com/css?family=Fira_sans');\n",
       "\n",
       "    div.cell{\n",
       "        width: 900px;\n",
       "        margin-left: 0% !important;\n",
       "        margin-right: auto;\n",
       "    }\n",
       "    div.text_cell code {\n",
       "        background: transparent;\n",
       "        color: #000000;\n",
       "        font-weight: 600;\n",
       "        font-size: 11pt;\n",
       "        font-style: bold;\n",
       "        font-family:  'Source Code Pro', Consolas, monocco, monospace;\n",
       "   }\n",
       "    h1 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "\t}\n",
       "\t\n",
       "    div.input_area {\n",
       "        background: #F6F6F9;\n",
       "        border: 1px solid #586e75;\n",
       "    }\n",
       "\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 30pt;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h2 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "        text-align: left;\n",
       "    }\n",
       "    .text_cell_render h2 {\n",
       "        font-weight: 200;\n",
       "        font-size: 16pt;\n",
       "        font-style: italic;\n",
       "        line-height: 100%;\n",
       "        color:#c76c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 1.5em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    } \n",
       "    h3 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h3 {\n",
       "        font-weight: 200;\n",
       "        font-size: 14pt;\n",
       "        line-height: 100%;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 2em;\n",
       "        display: block;\n",
       "        white-space: wrap;\n",
       "        text-align: left;\n",
       "    }\n",
       "    h4 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h4 {\n",
       "        font-weight: 100;\n",
       "        font-size: 14pt;\n",
       "        color:#d77c0c;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    h5 {\n",
       "        font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "    }\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 200;\n",
       "        font-style: normal;\n",
       "        color: #1d3b84;\n",
       "        font-size: 16pt;\n",
       "        margin-bottom: 0em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "        white-space: nowrap;\n",
       "    }\n",
       "    div.text_cell_render{\n",
       "        font-family: 'Fira sans', verdana,arial,sans-serif;\n",
       "        line-height: 125%;\n",
       "        font-size: 115%;\n",
       "        text-align:justify;\n",
       "        text-justify:inter-word;\n",
       "    }\n",
       "    div.output_subarea.output_text.output_pyout {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_subarea.output_stream.output_stdout.output_text {\n",
       "        overflow-x: auto;\n",
       "        overflow-y: scroll;\n",
       "        max-height: 50000px;\n",
       "    }\n",
       "    div.output_wrapper{\n",
       "        margin-top:0.2em;\n",
       "        margin-bottom:0.2em;\n",
       "}\n",
       "\n",
       "    code{\n",
       "      font-size: 70%;\n",
       "    }\n",
       "    .rendered_html code{\n",
       "    background-color: transparent;\n",
       "    }\n",
       "    ul{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li li{\n",
       "        padding-left: 0.2em; \n",
       "        margin-bottom: 0.2em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    ol{\n",
       "        margin: 2em;\n",
       "    }\n",
       "    ol li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.5em; \n",
       "    }\n",
       "    ul li{\n",
       "        padding-left: 0.5em; \n",
       "        margin-bottom: 0.5em; \n",
       "        margin-top: 0.2em; \n",
       "    }\n",
       "    a:link{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:visited{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:hover{\n",
       "       font-weight: bold;\n",
       "       color: #1d3b84;\n",
       "    }\n",
       "    a:focus{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    a:active{\n",
       "       font-weight: bold;\n",
       "       color:#447adb;\n",
       "    }\n",
       "    .rendered_html :link {\n",
       "       text-decoration: underline; \n",
       "    }\n",
       "    .rendered_html :hover {\n",
       "       text-decoration: none; \n",
       "    }\n",
       "    .rendered_html :visited {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :focus {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .rendered_html :active {\n",
       "      text-decoration: none;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "    } \n",
       "    hr {\n",
       "      color: #f3f3f3;\n",
       "      background-color: #f3f3f3;\n",
       "      height: 1px;\n",
       "    }\n",
       "    blockquote{\n",
       "      display:block;\n",
       "      background: #fcfcfc;\n",
       "      border-left: 5px solid #c76c0c;\n",
       "      font-family: 'Open sans',verdana,arial,sans-serif;\n",
       "      width:680px;\n",
       "      padding: 10px 10px 10px 10px;\n",
       "      text-align:justify;\n",
       "      text-justify:inter-word;\n",
       "      }\n",
       "      blockquote p {\n",
       "        margin-bottom: 0;\n",
       "        line-height: 125%;\n",
       "        font-size: 100%;\n",
       "      }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    scale:100,\n",
       "                        availableFonts: [],\n",
       "                        preferredFont:null,\n",
       "                        webFont: \"TeX\",\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# format the book\n",
    "%matplotlib inline\n",
    "import sys\n",
    "from __future__ import division, print_function\n",
    "import sys\n",
    "sys.path.insert(0,'../code')\n",
    "import book_format\n",
    "book_format.load_style('../code')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximate Bayesian Computation\n",
    "\n",
    "## The Variability Hypothesis\n",
    "\n",
    "I have a soft spot for crank science. Recently I visited Norumbega\n",
    "Tower, which is an enduring monument to the crackpot theories of Eben\n",
    "Norton Horsford, inventor of double-acting baking powder and fake\n",
    "history. But that’s not what this chapter is about.\n",
    "\n",
    "This chapter is about the Variability Hypothesis, which\n",
    "\n",
    "> “originated in the early nineteenth century with Johann Meckel, who\n",
    "> argued that males have a greater range of ability than females,\n",
    "> especially in intelligence. In other words, he believed that most\n",
    "> geniuses and most mentally retarded people are men. Because he\n",
    "> considered males to be the ’superior animal,’ Meckel concluded that\n",
    "> females’ lack of variation was a sign of inferiority.”\n",
    ">\n",
    "> From <http://en.wikipedia.org/wiki/Variability_hypothesis>.\n",
    "\n",
    "I particularly like that last part, because I suspect that if it turns\n",
    "out that women are actually more variable, Meckel would take that as a\n",
    "sign of inferiority, too. Anyway, you will not be surprised to hear that\n",
    "the evidence for the Variability Hypothesis is weak.\n",
    "\n",
    "Nevertheless, it came up in my class recently when we looked at data\n",
    "from the CDC’s Behavioral Risk Factor Surveillance System (BRFSS),\n",
    "specifically the self-reported heights of adult American men and women.\n",
    "The dataset includes responses from 154407 men and 254722 women. Here’s\n",
    "what we found:\n",
    "\n",
    "-   The average height for men is 178 cm; the average height for women\n",
    "    is 163 cm. So men are taller, on average. No surprise there.\n",
    "\n",
    "-   For men the standard deviation is 7.7 cm; for women it is 7.3 cm. So\n",
    "    in absolute terms, men’s heights are more variable.\n",
    "\n",
    "-   But to compare variability between groups, it is more meaningful to\n",
    "    use the coefficient of variation (CV), which is the standard\n",
    "    deviation divided by the mean. It is a dimensionless measure of\n",
    "    variability relative to scale. For men CV is 0.0433; for women it is\n",
    "    0.0444.\n",
    "\n",
    "That’s very close, so we could conclude that this dataset provides weak\n",
    "evidence against the Variability Hypothesis. But we can use Bayesian\n",
    "methods to make that conclusion more precise. And answering this\n",
    "question gives me a chance to demonstrate some techniques for working\n",
    "with large datasets.\n",
    "\n",
    "I will proceed in a few steps:\n",
    "\n",
    "1.  We’ll start with the simplest implementation, but it only works for\n",
    "    datasets smaller than 1000 values.\n",
    "\n",
    "2.  By computing probabilities under a log transform, we can scale up to\n",
    "    the full size of the dataset, but the computation gets slow.\n",
    "\n",
    "3.  Finally, we speed things up substantially with Approximate Bayesian\n",
    "    Computation, also known as ABC.\n",
    "\n",
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/variability.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean and standard deviation\n",
    "\n",
    "In Chapter [paintball] we estimated two parameters simultaneously using\n",
    "a joint distribution. In this chapter we use the same method to estimate\n",
    "the parameters of a Gaussian distribution: the mean, `mu`,\n",
    "and the standard deviation, `sigma`.\n",
    "\n",
    "For this problem, I define a Suite called `Height` that\n",
    "represents a map from each `mu, sigma` pair to its\n",
    "probability:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import thinkbayes\n",
    "\n",
    "class Height(thinkbayes.Suite, thinkbayes.Joint):\n",
    "\n",
    "    def __init__(self, mus, sigmas):\n",
    "        thinkbayes.Suite.__init__(self)\n",
    "\n",
    "        pairs = [(mu, sigma) \n",
    "                 for mu in mus\n",
    "                 for sigma in sigmas]\n",
    "\n",
    "        thinkbayes.Suite.__init__(self, pairs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`mus` is a sequence of possible values for `mu`;\n",
    "`sigmas` is a sequence of values for `sigma`. The\n",
    "prior distribution is uniform over all `mu, sigma` pairs.\n",
    "\n",
    "The likelihood function is easy. Given hypothetical values of\n",
    "`mu` and `sigma`, we compute the likelihood of a\n",
    "particular value, `x`. That’s what\n",
    "`EvalGaussianPdf` does, so all we have to do is use it:\n",
    "\n",
    "```python\n",
    "    # class Height\n",
    "\n",
    "        def Likelihood(self, data, hypo):\n",
    "            x = data\n",
    "            mu, sigma = hypo\n",
    "            like = thinkbayes.EvalGaussianPdf(x, mu, sigma)\n",
    "            return like\n",
    "```\n",
    "\n",
    "If you have studied statistics from a mathematical perspective, you know\n",
    "that when you evaluate a PDF, you get a probability density. In order to\n",
    "get a probability, you have to integrate probability densities over some\n",
    "range.\n",
    "\n",
    "But for our purposes, we don’t need a probability; we just need\n",
    "something proportional to the probability we want. A probability density\n",
    "does that job nicely.\n",
    "\n",
    "The hardest part of this problem turns out to be choosing appropriate\n",
    "ranges for `mus` and `sigmas`. If the range is too\n",
    "small, we omit some possibilities with non-negligible probability and\n",
    "get the wrong answer. If the range is too big, we get the right answer,\n",
    "but waste computational power.\n",
    "\n",
    "So this is an opportunity to use classical estimation to make Bayesian\n",
    "techniques more efficient. Specifically, we can use classical estimators\n",
    "to find a likely location for `mu` and `sigma`,\n",
    "and use the standard errors of those estimates to choose a likely\n",
    "spread.\n",
    "\n",
    "If the true parameters of the distribution are $\\mu$ and $\\sigma$, and\n",
    "we take a sample of $n$ values, an estimator of $\\mu$ is the sample\n",
    "mean, `m`.\n",
    "\n",
    "And an estimator of $\\sigma$ is the sample standard variance,\n",
    "`s`.\n",
    "\n",
    "The standard error of the estimated $\\mu$ is $s / \\sqrt{n}$ and the\n",
    "standard error of the estimated $\\sigma$ is $s / \\sqrt{2 (n-1)}$.\n",
    "\n",
    "Here’s the code to compute all that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def FindPriorRanges(xs, num_points, num_stderrs=3.0):\n",
    "\n",
    "    # compute m and s\n",
    "    n = len(xs)\n",
    "    m = numpy.mean(xs)\n",
    "    s = numpy.std(xs)\n",
    "\n",
    "    # compute ranges for m and s\n",
    "    stderr_m = s / math.sqrt(n)\n",
    "    mus = MakeRange(m, stderr_m, num_stderrs)\n",
    "\n",
    "    stderr_s = s / math.sqrt(2 * (n-1))\n",
    "    sigmas = MakeRange(s, stderr_s, num_stderrs)\n",
    "\n",
    "    return mus, sigmas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xs` is the dataset. `num_points` is the desired number of\n",
    "values in the range. `num_stderrs` is the width of the range on each\n",
    "side of the estimate, in number of standard errors.\n",
    "\n",
    "The return value is a pair of sequences, `mus` and\n",
    "`sigmas`.\n",
    "\n",
    "Here’s `MakeRange`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MakeRange(estimate, stderr, num_stderrs):\n",
    "    spread = stderr * num_stderrs\n",
    "    array = numpy.linspace(estimate-spread,\n",
    "                           estimate+spread,\n",
    "                           num_points)\n",
    "    return array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy.linspace` makes an array of equally spaced elements\n",
    "between `estimate-spread` and `estimate+spread`,\n",
    "including both."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Update\n",
    "\n",
    "Finally here’s the code to make and update the suite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'xs' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-16-b325b3abc99f>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mmus\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigmas\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mFindPriorRanges\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnum_points\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0msuite\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mHeight\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmus\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msigmas\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mUpdateSet\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMaximumLikelihood\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'xs' is not defined"
     ]
    }
   ],
   "source": [
    "mus, sigmas = FindPriorRanges(xs, num_points)\n",
    "suite = Height(mus, sigmas)\n",
    "suite.UpdateSet(xs)\n",
    "print(suite.MaximumLikelihood())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This process might seem bogus, because we use the data to choose the\n",
    "range of the prior distribution, and then use the data again to do the\n",
    "update. In general, using the same data twice is, in fact, bogus.\n",
    "\n",
    "But in this case it is ok. Really. We use the data to choose the range\n",
    "for the prior, but only to avoid computing a lot of probabilities that\n",
    "would have been very small anyway. With `num_stderrs=4`, the range is\n",
    "big enough to cover all values with non-negligible likelihood. After\n",
    "that, making it bigger has no effect on the results.\n",
    "\n",
    "In effect, the prior is uniform over all values of `mu` and\n",
    "`sigma`, but for computational efficiency we ignore all the\n",
    "values that don’t matter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The posterior distribution of CV\n",
    "\n",
    "Once we have the posterior joint distribution of `mu` and\n",
    "`sigma`, we can compute the distribution of CV for men and\n",
    "women, and then the probability that one exceeds the other.\n",
    "\n",
    "To compute the distribution of CV, we enumerate pairs of `mu`\n",
    "and `sigma`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def CoefVariation(suite):\n",
    "    pmf = thinkbayes.Pmf()\n",
    "    for (mu, sigma), p in suite.Items():\n",
    "        pmf.Incr(sigma/mu, p)\n",
    "    return pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use `thinkbayes.PmfProbGreater` to compute the probability that\n",
    "men are more variable.\n",
    "\n",
    "The analysis itself is simple, but there are two more issues we have to\n",
    "deal with:\n",
    "\n",
    "1.  As the size of the dataset increases, we run into a series of\n",
    "    computational problems due to the limitations of floating-point\n",
    "    arithmetic.\n",
    "\n",
    "2.  The dataset contains a number of extreme values that are almost\n",
    "    certainly errors. We will need to make the estimation process robust\n",
    "    in the presence of these outliers.\n",
    "\n",
    "The following sections explain these problems and their solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Underflow\n",
    "\n",
    "If we select the first 100 values from the BRFSS dataset and run the\n",
    "analysis I just described, it runs without errors and we get posterior\n",
    "distributions that look reasonable.\n",
    "\n",
    "If we select the first 1000 values and run the program again, we get an\n",
    "error in `Pmf.Normalize`:\n",
    "\n",
    "    ValueError: total probability is zero.\n",
    "\n",
    "The problem is that we are using probability densities to compute\n",
    "likelihoods, and densities from continuous distributions tend to be\n",
    "small. And if you take 1000 small values and multiply them together, the\n",
    "result is very small. In this case it is so small it can’t be\n",
    "represented by a floating-point number, so it gets rounded down to zero,\n",
    "which is called `**underflow**`. And if all probabilities in\n",
    "the distribution are 0, it’s not a distribution any more.\n",
    "\n",
    "A possible solution is to renormalize the Pmf after each update, or\n",
    "after each batch of 100. That would work, but it would be slow.\n",
    "\n",
    "A better alternative is to compute likelihoods under a log transform.\n",
    "That way, instead of multiplying small values, we can add up log\n",
    "likelihoods. `Pmf` provides methods `Log`,\n",
    "`LogUpdateSet` and `Exp` to make this process\n",
    "easy.\n",
    "\n",
    "`Log` computes the log of the probabilities in a Pmf:\n",
    "\n",
    "```python\n",
    "# class Pmf\n",
    "    def Log(self):\n",
    "        m = self.MaxLike()\n",
    "        for x, p in self.d.iteritems():\n",
    "            if p:\n",
    "                self.Set(x, math.log(p/m))\n",
    "            else:\n",
    "                self.Remove(x)\n",
    "```\n",
    "\n",
    "Before applying the log transform `Log` uses\n",
    "`MaxLike` to find `m`, the highest probability in\n",
    "the Pmf. It divide all probabilities by `m`, so the highest\n",
    "probability gets normalized to 1, which yields a log of 0. The other log\n",
    "probabilities are all negative. If there are any values in the Pmf with\n",
    "probability 0, they are removed.\n",
    "\n",
    "While the Pmf is under a log transform, we can’t use\n",
    "`Update`, `UpdateSet`, or `Normalize`.\n",
    "The result would be nonsensical; if you try, Pmf raises an exception.\n",
    "Instead, we have to use `LogUpdate` and\n",
    "`LogUpdateSet`.\n",
    "\n",
    "Here’s the implementation of `LogUpdateSet`:\n",
    "\n",
    "```python\n",
    "# class Suite\n",
    "    def LogUpdateSet(self, dataset):\n",
    "        for data in dataset:\n",
    "            self.LogUpdate(data)\n",
    "```\n",
    "\n",
    "`LogUpdateSet` loops through the data and calls\n",
    "`LogUpdate`:\n",
    "\n",
    "```python\n",
    "# class Suite\n",
    "    def LogUpdate(self, data):\n",
    "        for hypo in self.Values():\n",
    "            like = self.LogLikelihood(data, hypo)\n",
    "            self.Incr(hypo, like)\n",
    "```\n",
    "\n",
    "`LogUpdate` is just like `Update` except that it\n",
    "calls `LogLikelihood` instead of `Likelihood`, and\n",
    "`Incr` instead of `Mult`.\n",
    "\n",
    "Using log-likelihoods avoids the problem with underflow, but while the\n",
    "Pmf is under the log transform, there’s not much we can do with it. We\n",
    "have to use `Exp` to invert the transform:\n",
    "\n",
    "```python\n",
    "# class Pmf\n",
    "    def Exp(self):\n",
    "        m = self.MaxLike()\n",
    "        for x, p in self.d.iteritems():\n",
    "            self.Set(x, math.exp(p-m))\n",
    "```\n",
    "\n",
    "If the log-likelihoods are large negative numbers, the resulting\n",
    "likelihoods might underflow. So `Exp` finds the maximum\n",
    "log-likelihood, `m`, and shifts all the likelihoods up by\n",
    "`m`. The resulting distribution has a maximum likelihood of\n",
    "1. This process inverts the log transform with minimal loss of\n",
    "precision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log-likelihood\n",
    "\n",
    "Now all we need is `LogLikelihood`.\n",
    "\n",
    "```python\n",
    "# class Height\n",
    "\n",
    "    def LogLikelihood(self, data, hypo):\n",
    "        x = data\n",
    "        mu, sigma = hypo\n",
    "        loglike = scipy.stats.norm.logpdf(x, mu, sigma)\n",
    "        return loglike\n",
    "```\n",
    "\n",
    "`norm.logpdf` computes the log-likelihood of the Gaussian\n",
    "PDF.\n",
    "\n",
    "Here’s what the whole update process looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'suite' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-18-9dfe9138a6f7>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mLog\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mLogUpdateSet\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mExp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0msuite\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mNormalize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'suite' is not defined"
     ]
    }
   ],
   "source": [
    "suite.Log()\n",
    "suite.LogUpdateSet(xs)\n",
    "suite.Exp()\n",
    "suite.Normalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To review, `Log` puts the suite under a log transform.\n",
    "`LogUpdateSet` calls `LogUpdate`, which calls\n",
    "`LogLikelihood`. `LogUpdate` uses\n",
    "`Pmf.Incr`, because adding a log-likelihood is the same as\n",
    "multiplying by a likelihood.\n",
    "\n",
    "After the update, the log-likelihoods are large negative numbers, so\n",
    "`Exp` shifts them up before inverting the transform, which is\n",
    "how we avoid underflow.\n",
    "\n",
    "Once the suite is transformed back, the probabilities are “linear”\n",
    "again, which means “not logarithmic”, so we can use\n",
    "`Normalize` again.\n",
    "\n",
    "Using this algorithm, we can process the entire dataset without\n",
    "underflow, but it is still slow. On my computer it might take an hour.\n",
    "We can do better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A little optimization\n",
    "\n",
    "This section uses math and computational optimization to speed things up\n",
    "by a factor of 100. But the following section presents an algorithm that\n",
    "is even faster. So if you want to get right to the good stuff, feel free\n",
    "to skip this section.\n",
    "\n",
    "`Suite.LogUpdateSet` calls `LogUpdate` once for\n",
    "each data point. We can speed it up by computing the log-likelihood of\n",
    "the entire dataset at once.\n",
    "\n",
    "We’ll start with the Gaussian PDF:\n",
    "\n",
    "$$\\frac{1}{\\sigma \\sqrt{2 \\pi}} \\exp \\left[ -\\frac{1}{2} \\left( \\frac{x-\\mu}{\\sigma} \\right)^2 \\right]$$\n",
    "\n",
    "and compute the log (dropping the constant term):\n",
    "\n",
    "$$-\\log \\sigma -\\frac{1}{2} \\left( \\frac{x-\\mu}{\\sigma} \\right)^2$$\n",
    "\n",
    "Given a sequence of values, $x_i$, the total log-likelihood is\n",
    "\n",
    "$$\\sum_i -\\log \\sigma - \\frac{1}{2} \\left( \\frac{x_i-\\mu}{\\sigma} \\right)^2$$\n",
    "\n",
    "Pulling out the terms that don’t depend on $i$, we get\n",
    "\n",
    "$$-n \\log \\sigma - \\frac{1}{2 \\sigma^2} \\sum_i (x_i - \\mu)^2$$\n",
    "\n",
    "which we\n",
    "can translate into Python:\n",
    "\n",
    "```python\n",
    "# class Height\n",
    "    def LogUpdateSetFast(self, data):\n",
    "        xs = tuple(data)\n",
    "        n = len(xs)\n",
    "\n",
    "        for hypo in self.Values():\n",
    "            mu, sigma = hypo\n",
    "            total = Summation(xs, mu)\n",
    "            loglike = -n * math.log(sigma) - total / 2 / sigma**2\n",
    "            self.Incr(hypo, loglike)\n",
    "```\n",
    "\n",
    "By itself, this would be a small improvement, but it creates an\n",
    "opportunity for a bigger one. Notice that the summation only depends on\n",
    "`mu`, not `sigma`, so we only have to compute it\n",
    "once for each value of `mu`.\n",
    "\n",
    "To avoid recomputing, I factor out a function that computes the\n",
    "summation, and `**memoize**` it so it stores previously\n",
    "computed results in a dictionary (see\n",
    "<http://en.wikipedia.org/wiki/Memoization>):\n",
    "\n",
    "```python\n",
    "def Summation(xs, mu, cache={}):\n",
    "    try:\n",
    "        return cache[xs, mu]\n",
    "    except KeyError:\n",
    "        ds = [(x-mu)**2 for x in xs]\n",
    "        total = sum(ds)\n",
    "        cache[xs, mu] = total\n",
    "        return total\n",
    "```\n",
    "\n",
    "`cache` stores previously computed sums. The `try`\n",
    "statement returns a result from the cache if possible; otherwise it\n",
    "computes the summation, then caches and returns the result.\n",
    "\n",
    "The only catch is that we can’t use a list as a key in the cache,\n",
    "because it is not a hashable type. That’s why\n",
    "`LogUpdateSetFast` converts the dataset to a tuple.\n",
    "\n",
    "This optimization speeds up the computation by about a factor of 100,\n",
    "processing the entire dataset (154 407 men and 254 722 women) in less\n",
    "than a minute on my not-very-fast computer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ABC\n",
    "\n",
    "But maybe you don’t have that kind of time. In that case, Approximate\n",
    "Bayesian Computation (ABC) might be the way to go. The motivation behind\n",
    "ABC is that the likelihood of any particular dataset is:\n",
    "\n",
    "1.  Very small, especially for large datasets, which is why we had to\n",
    "    use the log transform,\n",
    "\n",
    "2.  Expensive to compute, which is why we had to do so much\n",
    "    optimization, and\n",
    "\n",
    "3.  Not really what we want anyway.\n",
    "\n",
    "We don’t really care about the likelihood of seeing the exact dataset we\n",
    "saw. Especially for continuous variables, we care about the likelihood\n",
    "of seeing any dataset like the one we saw.\n",
    "\n",
    "For example, in the Euro problem, we don’t care about the order of the\n",
    "coin flips, only the total number of heads and tails. And in the\n",
    "locomotive problem, we don’t care about which particular trains were\n",
    "seen, only the number of trains and the maximum of the serial numbers.\n",
    "\n",
    "Similarly, in the BRFSS sample, we don’t really want to know the\n",
    "probability of seeing one particular set of values (especially since\n",
    "there are hundreds of thousands of them). It is more relevant to ask,\n",
    "“If we sample 100,000 people from a population with hypothetical values\n",
    "of $\\mu$ and $\\sigma$, what would be the chance of collecting a sample\n",
    "with the observed mean and variance?”\n",
    "\n",
    "For samples from a Gaussian distribution, we can answer this question\n",
    "efficiently because we can find the distribution of the sample\n",
    "statistics analytically. In fact, we already did it when we computed the\n",
    "range of the prior.\n",
    "\n",
    "If you draw $n$ values from a Gaussian distribution with parameters\n",
    "$\\mu$ and $\\sigma$, and compute the sample mean, $m$, the distribution\n",
    "of $m$ is Gaussian with parameters $\\mu$ and $\\sigma / \\sqrt{n}$.\n",
    "\n",
    "Similarly, the distribution of the sample standard deviation, $s$, is\n",
    "Gaussian with parameters $\\sigma$ and $\\sigma / \\sqrt{2 (n-1)}$.\n",
    "\n",
    "We can use these sample distributions to compute the likelihood of the\n",
    "sample statistics, $m$ and $s$, given hypothetical values for $\\mu$ and\n",
    "$\\sigma$. Here’s a new version of `LogUpdateSet` that does it:\n",
    "\n",
    "```python\n",
    "    def LogUpdateSetABC(self, data):\n",
    "        xs = data\n",
    "        n = len(xs)\n",
    "\n",
    "        # compute sample statistics\n",
    "        m = numpy.mean(xs)\n",
    "        s = numpy.std(xs)\n",
    "\n",
    "        for hypo in sorted(self.Values()):\n",
    "            mu, sigma = hypo\n",
    "\n",
    "            # compute log likelihood of m, given hypo\n",
    "            stderr_m = sigma / math.sqrt(n)\n",
    "            loglike = EvalGaussianLogPdf(m, mu, stderr_m)\n",
    "\n",
    "            #compute log likelihood of s, given hypo\n",
    "            stderr_s = sigma / math.sqrt(2 * (n-1))\n",
    "            loglike += EvalGaussianLogPdf(s, sigma, stderr_s)\n",
    "\n",
    "            self.Incr(hypo, loglike)\n",
    "```\n",
    "\n",
    "On my computer this function processes the entire dataset in about a\n",
    "second, and the result agrees with the exact result with about 5 digits\n",
    "of precision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Robust estimation\n",
    "\n",
    "![Contour plot of the posterior joint distribution of mean and standard\n",
    "deviation of height for men in the\n",
    "U.S.](figs/variability_posterior_male.pdf)\n",
    "\n",
    "[fig.variability1]\n",
    "\n",
    "![Contour plot of the posterior joint distribution of mean and standard\n",
    "deviation of height for women in the\n",
    "U.S.](figs/variability_posterior_female.pdf)\n",
    "\n",
    "[fig.variability2]\n",
    "\n",
    "We are almost ready to look at results, but we have one more problem to\n",
    "deal with. There are a number of outliers in this dataset that are\n",
    "almost certainly errors. For example, there are three adults with\n",
    "reported height of 61 cm, which would place them among the shortest\n",
    "living adults in the world. At the other end, there are four women with\n",
    "reported height 229 cm, just short of the tallest women in the world.\n",
    "\n",
    "It is not impossible that these values are correct, but it is unlikely,\n",
    "which makes it hard to know how to deal with them. And we have to get it\n",
    "right, because these extreme values have a disproportionate effect on\n",
    "the estimated variability.\n",
    "\n",
    "Because ABC is based on summary statistics, rather than the entire\n",
    "dataset, we can make it more robust by choosing summary statistics that\n",
    "are robust in the presence of outliers. For example, rather than use the\n",
    "sample mean and standard deviation, we could use the median and\n",
    "inter-quartile range (IQR), which is the difference between the 25th and\n",
    "75th percentiles.\n",
    "\n",
    "More generally, we could compute an inter-percentile range (IPR) that\n",
    "spans any given fraction of the distribution, `p`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MedianIPR(xs, p):\n",
    "    cdf = thinkbayes.MakeCdfFromList(xs)\n",
    "    median = cdf.Percentile(50)\n",
    "\n",
    "    alpha = (1-p) / 2\n",
    "    ipr = cdf.Value(1-alpha) - cdf.Value(alpha)\n",
    "    return median, ipr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xs` is a sequence of values. `p` is the desired\n",
    "range; for example, `p=0.5` yields the inter-quartile range.\n",
    "\n",
    "`MedianIPR` works by computing the CDF of `xs`,\n",
    "then extracting the median and the difference between two percentiles.\n",
    "\n",
    "We can convert from `ipr` to an estimate of\n",
    "`sigma` using the Gaussian CDF to compute the fraction of the\n",
    "distribution covered by a given number of standard deviations. For\n",
    "example, it is a well-known rule of thumb that 68% of a Gaussian\n",
    "distribution falls within one standard deviation of the mean, which\n",
    "leaves 16% in each tail. If we compute the range between the 16th and\n",
    "84th percentiles, we expect the result to be `2 \\* sigma`. So\n",
    "we can estimate `sigma` by computing the 68% IPR and dividing\n",
    "by 2.\n",
    "\n",
    "More generally we could use any number of `sigmas`.\n",
    "`MedianS` performs the more general version of this\n",
    "computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def MedianS(xs, num_sigmas):\n",
    "    half_p = thinkbayes.StandardGaussianCdf(num_sigmas) - 0.5\n",
    "\n",
    "    median, ipr = MedianIPR(xs, half_p * 2)\n",
    "    s = ipr / 2 / num_sigmas\n",
    "\n",
    "    return median, s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, `xs` is the sequence of values; `num_sigmas` is the\n",
    "number of standard deviations the results should be based on. The result\n",
    "is `median`, which estimates $\\mu$, and `s`, which\n",
    "estimates $\\sigma$.\n",
    "\n",
    "Finally, in `LogUpdateSetABC` we can replace the sample mean\n",
    "and standard deviation with `median` and `s`. And\n",
    "that pretty much does it.\n",
    "\n",
    "It might seem odd that we are using observed percentiles to estimate\n",
    "$\\mu$ and $\\sigma$, but it is an example of the flexibility of the\n",
    "Bayesian approach. In effect we are asking, “Given hypothetical values\n",
    "for $\\mu$ and $\\sigma$, and a sampling process that has some chance of\n",
    "introducing errors, what is the likelihood of generating a given set of\n",
    "sample statistics?”\n",
    "\n",
    "We are free to choose any sample statistics we like, up to a point:\n",
    "$\\mu$ and $\\sigma$ determine the location and spread of a distribution,\n",
    "so we need to choose statistics that capture those characteristics. For\n",
    "example, if we chose the 49th and 51st percentiles, we would get very\n",
    "little information about spread, so it would leave the estimate of\n",
    "$\\sigma$ relatively unconstrained by the data. All values of\n",
    "`sigma` would have nearly the same likelihood of producing\n",
    "the observed values, so the posterior distribution of `sigma`\n",
    "would look a lot like the prior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Who is more variable?\n",
    "\n",
    "![Posterior distributions of CV for men and women, based on robust\n",
    "estimators.](figs/variability_cv.pdf)\n",
    "\n",
    "[fig.variability3]\n",
    "\n",
    "Finally we are ready to answer the question we started with: is the\n",
    "coefficient of variation greater for men than for women?\n",
    "\n",
    "Using ABC based on the median and IPR with `num_sigmas=1`, I computed\n",
    "posterior joint distributions for `mu` and\n",
    "`sigma`. Figures [fig.variability1] and  [fig.variability2]\n",
    "show the results as a contour plot with `mu` on the x-axis,\n",
    "`sigma` on the y-axis, and probability on the z-axis.\n",
    "\n",
    "For each joint distribution, I computed the posterior distribution of\n",
    "CV. Figure [fig.variability3] shows these distributions for men and\n",
    "women. The mean for men is 0.0410; for women it is 0.0429. Since there\n",
    "is no overlap between the distributions, we conclude with near certainty\n",
    "that women are more variable in height than men.\n",
    "\n",
    "So is that the end of the Variability Hypothesis? Sadly, no. It turns\n",
    "out that this result depends on the choice of the inter-percentile\n",
    "range. With `num_sigmas=1`, we conclude that women are more variable,\n",
    "but with `num_sigmas=2` we conclude with equal confidence that men are\n",
    "more variable.\n",
    "\n",
    "The reason for the difference is that there are more men of short\n",
    "stature, and their distance from the mean is greater.\n",
    "\n",
    "So our evaluation of the Variability Hypothesis depends on the\n",
    "interpretation of “variability.” With `num_sigmas=1` we focus on people\n",
    "near the mean. As we increase `num_sigmas`, we give more weight to the\n",
    "extremes.\n",
    "\n",
    "To decide which emphasis is appropriate, we would need a more precise\n",
    "statement of the hypothesis. As it is, the Variability Hypothesis may be\n",
    "too vague to evaluate.\n",
    "\n",
    "Nevertheless, it helped me demonstrate several new ideas and, I hope you\n",
    "agree, it makes an interesting example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "There are two ways you might think of ABC. One interpretation is that it\n",
    "is, as the name suggests, an approximation that is faster to compute\n",
    "than the exact value.\n",
    "\n",
    "But remember that Bayesian analysis is always based on modeling\n",
    "decisions, which implies that there is no “exact” solution. For any\n",
    "interesting physical system there are many possible models, and each\n",
    "model yields different results. To interpret the results, we have to\n",
    "evaluate the models.\n",
    "\n",
    "So another interpretation of ABC is that it represents an alternative\n",
    "model of the likelihood. When we compute `$\\mathrm{p}(D|H)$`,\n",
    "we are asking “What is the likelihood of the data under a given\n",
    "hypothesis?”\n",
    "\n",
    "For large datasets, the likelihood of the data is very small, which is a\n",
    "hint that we might not be asking the right question. What we really want\n",
    "to know is the likelihood of any outcome like the data, where the\n",
    "definition of “like” is yet another modeling decision.\n",
    "\n",
    "The underlying idea of ABC is that two datasets are alike if they yield\n",
    "the same summary statistics. But in some cases, like the example in this\n",
    "chapter, it is not obvious which summary statistics to choose.\n",
    "\n",
    "You can download the code in this chapter from\n",
    "<http://thinkbayes.com/variability.py>. For more information see\n",
    "Section [download]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "An “effect size” is a statistic intended to measure the difference\n",
    "between two groups (see <http://en.wikipedia.org/wiki/Effect_size>).\n",
    "\n",
    "For example, we could use data from the BRFSS to estimate the difference\n",
    "in height between men and women. By sampling values from the posterior\n",
    "distributions of $\\mu$ and $\\sigma$, we could generate the posterior\n",
    "distribution of this difference.\n",
    "\n",
    "But it might be better to use a dimensionless measure of effect size,\n",
    "rather than a difference measured in cm. One option is to use divide\n",
    "through by the standard deviation (similar to what we did with the\n",
    "coefficient of variation).\n",
    "\n",
    "If the parameters for Group 1 are $(\\mu_1, \\sigma_1)$, and the\n",
    "parameters for Group 2 are $(\\mu_2, \\sigma_2)$, the dimensionless effect\n",
    "size is\n",
    "$$\\frac{\\mu_1 - \\mu_2}{(\\sigma_1 + \\sigma_2)/2}$$ \n",
    "\n",
    "Write a function that takes joint distributions of `mu` and\n",
    "`sigma` for two groups and returns the posterior distribution\n",
    "of effect size.\n",
    "\n",
    "Hint: if enumerating all pairs from the two distributions takes too\n",
    "long, consider random sampling."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
